library(SemiParBIVProbit)
library(Matching)
sim.func <- function(n.rep=100, n=500, d.unob = 1, beta = 0){

match.res <- matrix(NA, n.rep, 3)
true <- rep(NA, n.rep)


for(i in 1:n.rep){

## measured potential confounders
z1 <- runif(n)

## unmeasured confounder
if(d.unob==1){

z2 <- rnorm(n)       
                     
}                     
                     
if(d.unob==2){

z2 <- rt(n, df=4)       
} 

if(d.unob==3){

z2 <- rchisq(n, df=1)       
} 

if(d.unob==4){

z2 <- runif(n, -3,3)       
} 

  

z4 <- runif(n)                     
z3 <- rbinom(n,1,0.5) # IV


## non-linearities
f1   <- function(x) cos(pi*2*x) + sin(pi*x)
f2   <- function(x) x + exp(-30*(x-0.5)^2) 

## treatment assignment

prob.treated <-	plogis(-0.5 + f1(z1) + beta*z2 + 3*z3 + 1.3*z4) 
x <- rbinom(n, 1, prob.treated)

## potential outcomes
p0 <- plogis(-3.5 + f2(z1) + 2*z2 -0.8*z4)                    
p1 <- plogis( 0.5 + f2(z1) + 2*z2 -0.8*z4)                   


y0 <- rbinom(n, 1, p0)
y1 <- rbinom(n, 1, p1)

## observed outcomes
y <- y0
y[x==1] <- y1[x==1] 	     	

table(y)
table(x)
true[i] <- mean(y1)-mean(y0)

## put everything in a data frame
examp.data <- data.frame(z1, z2, z3, z4, x, y)


# Estimate propensity score and save the linear predictor of 
# treatment assignment 

tr.eq <- x ~ s(z1) + z3 + s(z4)

pscore_l <- gam(tr.eq, data = examp.data, binomial(link = "logit"))
pscore_p <- gam(tr.eq, data = examp.data, binomial(link = "probit"))
pscore_c <- gam(tr.eq, data = examp.data, binomial(link = "cloglog"))

pscore_est_l <- predict(pscore_l)

pscore_est_p <- predict(pscore_p)

pscore_est_c <- predict(pscore_c)

# The covariates we want to match on
X_l = cbind(pscore_est_l, z1, z4)
X_p = cbind(pscore_est_p, z1, z4)
X_c = cbind(pscore_est_c, z1, z4)

#The covariates we want to obtain balance on
BalanceMat_l <- X_l
BalanceMat_p <- X_p
BalanceMat_c <- X_c

# Let's call GenMatch() to find the optimal weight to give each
# covariate in 'X' so as we have achieved balance on the #covariates in
# 'BalanceMat'

genout_l <- GenMatch(Tr=x, X=X_l, BalanceMatrix=BalanceMat_l, estimand="ATE", M=1)

genout_p <- GenMatch(Tr=x, X=X_p, BalanceMatrix=BalanceMat_p, estimand="ATE", M=1)

genout_c <- GenMatch(Tr=x, X=X_c, BalanceMatrix=BalanceMat_c, estimand="ATE", M=1)

# Now that GenMatch() has found the optimal weights, let's #estimate
# our causal effect of interest using those weights
#

mout_l <- Match(Y=y, Tr=x, X=X_l, estimand="ATE", Weight.matrix=genout_l)
mout_p <- Match(Y=y, Tr=x, X=X_p, estimand="ATE", Weight.matrix=genout_p)
mout_c <- Match(Y=y, Tr=x, X=X_c, estimand="ATE", Weight.matrix=genout_c)



match.res[i,] <- c(mout_l$est, mout_p$est, mout_c$est)
                   
                 
                     
                     
                     
print(i)

}

   

true.v <- mean(true)
bias <- abs((colMeans(match.res[,1:3], na.rm = T)-true.v)/true.v)*100
names(bias) <- c(         "matchingL",
                          "matchingP",
                          "matchingC"
 )
			  
			  
			  
mse <- apply((match.res[,1:3] - true.v)^2, 2, mean, na.rm=TRUE)
names(mse) <- c  ( "matchingL",
                          "matchingP",
                          "matchingC"
			  )


results <- list(bias=bias, mse= mse, match.res=match.res, true.v= true.v)

results

}


set.seed(0)
results <- sim.func(n.rep = 250, n = 1000, d.unob = 2, beta = -0.85)

save.image("resbetaminus0.85unob2.RData")



bias <-results$bias
round(bias,2)

rmse <-sqrt(results$mse)
round(rmse,2)


